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BLACK HOLE-NEUTRON STAR MERGERS WITH A HOT NUCLEAR EQUATION OF STATE: 
OUTFLOW AND NEUTRINO-COOLED DISK FOR A LOW-MASS, HIGH-SPIN CASE 
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ABSTRACT 

Neutrino emission significantly affects the evolution of the accretion tori formed in black hole-neutron star 
mergers. It removes energy from the disk, alters its composition, and provides a potential power source for a 
gamma-ray burst. To study these effects, simulations in general relativity with a hot microphysical equation of 
state and neutrino feedback are needed. We present the first such simulation, using a neutrino leakage scheme 
for cooling to capture the most essential effects and considering a moderate mass ( 1 .4 Mq neutron star, 5.6 Mq 
black hole), high spin (black hole J/M^ = 0.9) system with the Kq = 220 MeV Lattimer-Swesty equation of 
state. We find that about 0.08 Mq of nuclear matter is ejected from the system, while another 0.3 Mq forms 
a hot, compact accretion disk. The primary effects of the escaping neutrinos are (/) to make the disk much 
denser and more compact, (//) to cause the average electron fraction Yg of the disk to rise to about 0.2 and 
then gradually decrease again, and (///) to gradually cool the disk. The disk is initially hot (T ^ 6 MeV) and 
luminous in neutrinos {L^, ^ 10^"* erg s"'), but the neutrino luminosity decreases by an order of magnitude over 
50 ms of post-merger evolution. 

Subject headings: Black hole physics - Gamma-ray burst: general - Neutrinos - Stars: neutron 



\. INTRODUCTION 

Much of the interest in black hole-neutron star (BHNS) 
mergers arises from their potential to solve two important 
problems in contemporary astrophysics. First, it is possible 
that such events can produce short hard gamma-ray bursts. 
Second, they could significantly contribute to the abundances 
of r-process nuclei observed in the solar system. 

Short hard gamma-ray bursts emit characteristic luminosi- 
ties of L ^ 10^*'"^^/(47r)ergs"' steradian"', with typical dura- 
tions of ^1 s, and peak photon energies of hv ^ 1 MeV. Their 
spectra are nont hermal, and va ry in brightness on a timescale 
of Af < 10 ms ( Nakar|[2007| l. Combined observations of a 
nonthermal spectrum arid a short time variability can be ex- 
plained by an ultra-relativistic jet (F > 30) launche d from a 
disk experiencing rapid a ccretion onto a black hole ( Narayan 
|erari|T992| |Nakar|p07l ). BHNS mergers are plausible pro- 
genitors for such a system. 

BHNS mergers may also produce unbound neutron-rich 
ejecta, which would provide an r-process nucleosynthesis 
site ne eded to explain the observed heavy element abun- 
dances (|Lattimer & Schramm|1974[|Fr eiburghaus et al.|1999 



or from the shock that would form when the outflow hits a suf - 
ficiently dense interstellar medium ( [Metzger & Berger|2012[ ). 

BHNS mergers and their aftermath involve relativistic grav- 
ity, magnetohydrodynamics, nuclear physics, and neutrino 
radiation. In addition, a large space of binary parameters 
must be explored, since a wide range of premerger black hole 



mass and spin va lues are plausible (Belczynski et al. 2008 



Ozel et al. 



2010 1. Early attempts to explore this parameter 



space in general rel ativity have used polytropic (TanigucMl 
et al. 2005 ; Shibata & Uryu 2006; Etienne et al. 2008 200^ 
Duez et al. 2008 2010; Foucart et al. 2011, "' 



etal.|2012 
toku et al. 



,2012j |Etienne 



Foucart etal. 2013)orpiecewise-polytropic (Kyu- 
01 0||2011( [Lackey et al.|2013) equations of state 



(EOS) to describe the neutron star matter, with the thermal 
effects being modeled by a simple F-law. Because the mat- 
ter does not heat until after tidal disruption, use of these F- 
law EOSs may be adequate for the inspiral and very early 
tidal disruption phases, which produce nearly all of the grav- 
itational wave signal. From these parameter s pace s t udies, 
a gene ral picture has beg un to emerge (Faber '2009'; 'Duez 
2010. :^imarale et al.|20Tl ; Shibata & Taniguchi 201 1 ; Fou 



[Arnould et al. 2007 ; Korobkin et al. 20l2|. Additionally such 
outflows may be observable from the radioactive decay that 
results after unstable heavy isotopes are formed (a 'kilonova') 
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|cart||2012'] [ In the best-understood binary parameter space 
of spin {a* = (/bh/Mbh^)(c/G) = [-0.5,0.9]), NS compact- 
ness (C = (Mns/-Rns)(G/c2) - [0.13,0.20]), and mass ra- 
tio iq = Mbh/A^ns ^ [1,7]) — where G is Newton's gravita- 
tional constant and c is the speed of light — disks of significant 
masses appear to be formed in binaries with low mass ratios, 
low NS compactness, or high prograde black hole spins. 

After the merger, neutrinos will begin to play an impor- 
tant role in the evolution. Following disruption, the fluid is 
shocked to temperatures of ~10 MeV and begins to radiate 
copiou sly. In simulations of accre tion disk-black hole sys- 
tems (Shibata & Sekiguchi"2012 l, and binary neutron star 
mergers (Ruffert et al. 1997; Rosswog & Liebendorfer|2003| 



Sekiguchi et al.|[201 1 1, which produce disks similar to those 
of BHNS mergers, total neutrino luminosities are of order 
10^'"^^ ergs"'. BHNS simulations using F-law EOSs indicate 
that disks last at most several hundred milliseconds (Etienne 



et al. 2009; |Duez et al. 2010; jKyutoku et al.|[20T()) |Foucart| 
et al.„2011i . During this time, neutrmos may carry away a 
significant amount of energy from the accretion disk, causing 
it to cool and contract. Additionally, neutrinos, unlike pho- 
tons, change the composition of their source. In the present 
case, initially neutron-rich neutron star material releptonizes 
as it expands into a disk through an imbalance of electron- 
and positron-capture reactions. These composition changes 
can play a significant role in BHNS post-merger dynamics 
in a number of ways. Changes to the electron fraction alter 
the luminosity of each neutrino species. Lepton number gra- 
die nts affect pressure forces and can drive convection (|Lee| 
|et al. 2005 ). Phase transitions, e.g. recombination from pure 
nucleonic matter into alpha particles and heavy nuclei, can 
release or absorb energy as density and temperature change. 
A temperature- and composition-dependent EOS is needed to 
capture all of these effects. 

Neutrino emission may also play a crucial role in powering 
the relativistic outflow needed for a GRB. It could do this by 
depositing energy in the funnel region along the BH spin axis 
via neutrino pair annihilation {lyi? — > e~e'^). Depending on the 
shape of the emitting region, its proximity to the black hole, 
and its emission spectrum, the annihilation proces s could pro- 
ceed with an efficiency as high as 0.2%-0.5% (jRuffert et al. 
1997; Birk l et al.|[2007l [Dessart et al1|2009l [Harikae et al" 
2010p. This would indicate that neutrino pair annihilation 
may only be able to power low-energy GRBs. However, ac- 
curate relativistic simulations of the fluid and neutrino fields 
are needed to make reliable statements about the efficiency of 
the neutrino pair annihilation process in BHNS mergers. 

Because the disk spans optically thick and optically thin 
regimes, a full 6-dimensional evolution of the neutrino fields 
(1 energetic + 3 spatial + 2 angular dimensions) is needed to 
completely describe the coupling between radiation and mat- 
ter This is impossible with current resources. Fortunately, 
many of the essential features of the radiation can be captured 
using a simple leakage scheme. Rather than performing ac- 
tual radiation transport, leakage schemes remove energy and 
alter lepton number at rates based on the local free-emission 
and diffusion rates. Leakage certainly neglects some arguably 
important effects (e.g. neutrino absorption-driven winds), but 
it captures the basic energetics and composition drift of the 
post-merger system. (See e.g. Ott et al. 2012]for an analysis 
of leakage in core-collapse supernovae simulations.) 

In Newtonian gravity, microphysical EOSs and neutrino 
leakage s chemes have been used in simulations of binary neu- 
ti-on star CRuffert et al.|1996|[T997l|Rosswog & Davies 



Rosswog & Liebend6 rfer|2003| l and BHNS ( |Janka et al 



2002 



1999 



osswog et al.||2004) mergers. Leakage schemes have re 



cently been incorporated into three-dimensional general rela- 
tivistic simulations of stella r core collapse (Sekiguchi & Shi- 
bate 2011 Ott et al.]2012i and binary neu tron star mergers 
( [Sekiguchi et al.||2011[|Kiuchi et al.||2012| l. Simulations of 
BHNS mergers with a micr ophysical E OSs but no neutrino 
feedback were presented by |Duez et al. !(2010). The first rel- 
ativistic BHNS simulations with both a microphysical EOS 
and neutrino feedback are presented here. 

In this paper, we introduce a neutrino leakage scheme 



into the Spectral Einstein Coda^ evolving the coupled 
hydrodynamics-Einstein field equations with leakage, and ap- 
ply it to one BHNS case, a moderately low-mass (Mns = 
IAMq, Mbh = 5.6M0), high black hole spin (a* = 0.9) sys- 
tem. We follow six orbits of inspiral and evolve for 50ms 
past merger. We employ the finite-temperature EOS of Lat- 



timer & Swesty] ( 1991| l using a compressibility parameter 
Kq = 220 MeV that yields neutron stars consistent with ex- 



isting obser vations (see e.g. Demorest et al.|20T0 



and Steiner 
[eraiT|[20l0l — though recent work by |Uuillot et al.|[2013| pre- 
dicts smaller neutron star radii). 

We find an initially large postmerger disk of mass Mo ^ 
0.3 M0 which accretes rapidly for the first 30 ms and then set- 
tles to a slow mass depletion rate, from which we can esti- 
mate a disk lifetime of T^ep > 0.2 s. The neutrino luminos- 
ity peaks 10 ms after merger at ^10'''* ergs"' but drops to 
^2 X 10^^ ergs"' by the end of the 50 ms post-merger evo- 
lution, at which time it continu es to drop. Assu ming a conser- 
vative efficiency of 0.2% (Ruffert et al.||1997 ), we estimate 
the average work done by neutrino pair annihilation to be 
2i/p ^ 10^' ergs"', for a total energy available to drive a jet 
of E^fy ^ 10^'^ erg. We find a large unbound outflow from the 
system, of mass M^ ^ 0.08 M©, with mildly relativistic ve- 
locities, leading to an available kinetic energy that could be as 
high as Ef.j ^ 10^^ erg. 

The disk itself displays an interesting evolution during this 
time. A very high energy region — both hot and having super- 
Keplerian kinetic energies^-develops in the innermost disk. 
This configuration is unstable, and nonaxisymmetric pertur- 
bations persist in the inner disk for many dynamical times. 
Away from the edges, the disk does roughly settle to station- 
ary axisymmetry, and evolution is driven primarily by secular 
radiation and disk depletion effects. Neutrino cooling leads 
to a much higher-density, and somewhat lower-entropy, disk 
structure than seen in a comparison run without cooling. The 
average temperature decreases in the disk from a maximum of 
^6 MeV at 15 ms after merger to ^4 MeV at 50 ms. The av- 
erage specific entropy (~0. 1 ^b baryon"' during inspiral) stays 
near ^8 ^b baryon"' from 5 ms to the end of the evolution. 
The disk is sufficiently hot and dense to remain opaque to 
neutrinos of all species. The average electron fraction rises 
from 0.07 in the initial neutron star to around 0.2 at 20 ms 
after merger. At this time, electron neutrino and antineutrino 
number emission roughly balance. Afterwards, the electron 
fraction decreases gradually as the disk cools and the balance 
of ff and Dg emission adjusts. The disk continues to cool and 
flatten, so that by the end of our simulation it probably no 
longer has the energy to power a GRB via a neutrino mech- 
anism alone. Longer-lasting high-power energy release may 
come from physics not included in this simulation, particu- 
larly magnetic fields. 

This paper is organized as follows. In Section |2] we dis- 
cuss the numerical methods used to simulate the mergers; in 
subsections |2.1[[2!2l and |2.3| we describe our treatment of the 
nuclear EOS, initial data, and neutrino leakage, respectively. 
In Section I3] we give a summary of the evolution, examining 
several measures of convergence. In Section l4j we analyze 
outflows. In Section |5] we review the relevant timescales for 
the accretion disk. In Sections |6] and |7] we analyze the accre- 
tion disk in the epochs of formation and neutrino-driven evo- 



http : //www. black-holes . org/SpEC .html 



lution, respectively. We summarize our conclusions in Sec- 
tionE 

Throughout the rest of this paper, unless otherwise noted, 
we adopt geometric units in which G = c = 1, and we use 
Latin indices {j,k) to represent the three spatial coordinates, 
and Greek indices (a, f3) to represent the four spacetime coor- 
dinates. 

2. METHODS 

We evolve the coupled general relativistic-hydrodyamics 
system using the Spectral Einstein Code (SpEC). We em- 
ploy our standard two-grid pseudospectral/fini te difference 
approach, described in detail in earlier papers d Scheel et al. 
2006} |Duez et al.|2008[ [Hemberger eFal 2012; F oucart et aH 
2013') and briefly sumimarized here. The spacetime is de- 
scribed by the 4-metric and its time derivative, which we 
evolve in the generalized harmonic formulat ion using a mul - 



tidomain pseudospectral algorithm (Lindblom et al. 2006|l. 



The coordinates are evolved by enforcing the 'frozeir gauge 



condition described in Appendix A.l of Foucart et al. ( 2013) l. 
The fluid is described by the hydrodynamic fields, which 
we evolve in conservative form using shock-capturing finite- 
difference methods. The metric and fluid are evolved on sep- 
arate computational domains; fluid source terms for the met- 
ric evolution and metric source terms for the fluid evolution 
are acquired by interpolation between domains. The fluid do- 
main is a uniform Cartesian grid covering the non-vacuum 
upper hemisphere, with the lower hemisphere fluid quantities 
set by an assumed equatorial symmetry. The black hole sin- 
gularity is handled by excising a region inside the apparent 
horizon, i.e. by not placing colocation points there. The fluid 
grid does have points inside the excised region, but the fields 
at these points are not evolved, and one-sided stencils are used 
to evolve points next to the excision boundary. 

The fluid is described by the baryonic rest mass density 
p, temperature T, electron fraction (electrons per baryon) Y,,, 
and 3-velocity v'. The pressure P and specific internal energy 
e are computed from p, T, and Y^ using an EOS table (see 
Section [ZT] below). From these, one can compute the spe- 
cific enthalpy h= 1 + e + P/p. From v', one can compute the 
Lorentz factor W = au', a being the lapse and u' being the 
time component of the 4-velocity. The hydrodynamic evo- 
lution equations take the form of conservation laws for a set 
of 'conservative' variables. These include a density variable, 
p* = sJ^fN p, an energy variable, t = p^^ihW -\)- ^P , a mo- 
mentum variable 5, = p^hui, and a composition variable p^Ff , 
where 7 is the determinant of the 3-metric, and «,- are the co- 
variant components of the 4-velocity. Note that we do not 
include neutrino pressure or energy in the definitions of f and 
Si\ the neutrino radiation field appears only through source 
terms in the evolution equations (see Section 2.3 1. 



To carry out a timestep, we begin by computing the con- 
servative variabl e fluxes o n cell faces using 5t h-orde r WENO 
reconstruction (Liu et al.l 1994; Jiang & Shu||1996|l of p, T, 



and M , and an HLL approximate Riemann solver ( |A. Harten| 
|1983| l. From these cell boundary fluxes, the conservative vari- 
ables are advanced forward in time. From the evolved conser- 
vative variables, the 'primitive' variables p, T, and v' must 
then be recovered. This amounts to a two-dimensional root- 
finding problem to find the values of W and T that reproduce 
the evolved conservative variables. For some combinations 



of conservative variables, no corresponding set of primitive 
variables exists. For this reason, relativistic hydrod ynamics 
codes often impose the condition S^ < f(f + 2p*) (Etiennel 
iet al. 2008 ). This inequality assumes h > 1 for all physical 
values of p and T. One feature of nuclear-based EOSs is that 
the internal energy becomes negative at sufficiently low den- 
sities and temperatures. For an EOS with a specific enthalpy 
minimum of hmin, the invertibility condition becomes 



^ < '^max ■ 



■■f(f + 2p,) + (l-hlj 



(1) 



To avoid divisions by zero, we impose a density floor of 
6 X lO-'gcm"''. Gas near the surface, at points with densi- 
ties more than a few decades below the maximum, cannot 
be expected to be evolved accurately, so high temperatures 
tend to develop there. This is not a serious problem, but it 
can slow down the primitive variable recovery and introduce 
some anomalous neutrino luminosity. We therefore impose 
a maximum value on the temperature in the low-density re- 
gion outside the star, tail, and disk (lOMeV for densities less 
than lO"'* of the instantaneous maximum density). By running 
segments of the evolution with different density thresholds for 
temperature manipulation, we confirm that its main effect is to 
reduce noise in the neutrino luminosity during tidal disruption 
— noise that is dwarfed by the post-disruption signal anyway. 

2.1. Equation of State 

In this work we describe the fluid with the Lattimer & 
Swesty (hereafter LS) EOS. This EOS is derived from a 
compressible liquid-drop model with a Skyrme nuclear force 
and includes contributions from free nu cleons, alpha parti- 
cles, a nd a single type of heavy nucleus (Lattimer & Swesty 
1991|l. We set the nuclear incompressibility parameter, Kq 



to 220 MeV and the symmetry energy Sy to 29.3 MeV. Elec- 
trons, positrons, a nd photons are added using the routines of 
[Timmes & Arnett| ([T999). 

We employ the LS EOS in our evolutions as a table of the 
following thermodynamic and composition quantities: inter- 
nal energy; pressure; sound speed; neutron, proton, and elec- 
tron chemical potentials; and neutron, proton, alpha particle 
and characteristic heavy nucleus mass fractions (along with 
the characteristic heavy nucleus's average mass and charge 
number). We store the table with the following ranges and 
resolutions: pG [10^, 10"']gcm~-', N = 250, log spacing; 
T e [0.01,251] MeV, A^= 120, l og spacing; 7, € [0 035,0.53], 
A^ = 100, linear spacing. As in [O'Connor & Ott| ( 12010 ) we 
perform tri-linear interpolations for intermediate values. The 
original table and access routines are available at http : 
//www. stellarcollapse . org During evolution, T at 
a fluid grid point may evolve to a value below 0.01 MeV, or Yg 
may evolve below 0.035. In this case, the out-of-range vari- 
able is reset to the relevant table minimum. Continuous ex- 
trapolations of P and e outside the tabular bounds are defined 
for convenience in primitive variable recovery, but because T 
and Yg are immediately reset to remain within tabular bounds, 
these extrapolations do not influence the evolution. 

Just below nuclear saturation density a phase transition oc- 
curs between nuclei at low densities and uniform nuclear mat- 
ter at high densities. The LS EOS captures this transition 
smoothly in free energy, but other parameters suffer from 
jumps between the two phases. During the inversion from 
conservative to primitive variables, we solve for temperature 



an equation involving enthalpy and pressure. The inversion, 
therefore, is difficult when either h{T) or P(T) are not smooth. 
We find that the sharp transition is most disruptive to inver- 
sion at densities near p ^ 10'"* gem"-', and temperatures near 
r ^ 1 MeV. We resolve the issue by 'polishing' the table in the 
temperature dimension, applying a gaussian smoothing kernel 
of width 0.05 MeV to the stored values of internal energy and 
pressure. 

2.2. Initial data 

We choose initial parameters that we expect to yield a mas- 
sive accretion disk, and high accretion efficiencies, giving us 
an upper bound on energetics. The binary is characterized 
by a low mass ratio of q = 4, a high prograde aligned spin 
of a* = 0.9, and an initial orbital separation yielding 6 orbits 
of inspiral before disruption. In addition, we set the star's 
initial temperature to 0.01 MeV, and choose a Y,, profile that 
enforces neutrinoless /3-equilibrium. The star's baryonic rest 
mass is 1.55M0. In isolation, it would have a gravitational 
mass of 1.4OM0 and an areal radius of 12.7 km. 

To compute initial data we solve the extended conformal 
thin sandwich equations using multi-domain pseudospectral 
methods, as detailed in Foucart et al. ( 2008 ), using the Spec - 
tral Elliptic Solver ( Pfeiffer,2003nPfeiffer et al.|2003|[2005] l. 
The numerical domain is decomposed into a set of touch- 
ing but not overlapping 'subdomains', comprised of spheri- 
cal shells, filled spheres, cylinders, and distorted cubes, ar- 
ranged to reflect the symmetries in the configuration. We use 
the same method to re present the metric in evolutions (see 
e.g. Figures 1 and 2 in |Foucart et al.|[2013| l. One difficulty 
arising within a spectral framework is the spurious Gibbs os- 
cillations that can occur at discontinuities, for example at the 
surface of the star Thus, a critical technique in our method in- 
volves capturing the surface at a subdomain boundary, where 
the abutting spectral domains have no difficulty representing 
the no n-smooth field. Fo r the polytropes used in our previous 
work (|Duez et al.|[20T0l |Foucart et al.1[20n| [20T2l |Foucart| 
|2012[ [Foucart et al.||2013| ) the density profile is smooth and 
easy to resolve. However new difficulties present themselves 
with the complexity of the stellar structure derived from the 
LS EOS. 

Polytropes have density profiles that fall off as p ex (r-r^,)" 
close to the surface (Gundlach & Please 2009), where r^ is 
the stellar radius, and n is related to the polytropic index by 
r = l + l/n. Thus, density profiles of F = 2 polytropes behave 
linearly at the surface. To understand the effect of realistic 
EOSs on stellar surfaces we may use the effective adiabatic 
index, defined f = (dlogP/dlogp)s. The LS EOS, at low 
temperature and Y^ in neutrinoless /3-equilibrium, is extremely 
stiff near the central density of a 1.4Mq star, with T ~ 7/2. 
It begins to soften dramatically at ~' 10'"* gem"-', to F ^ 1/2. 
At even lower densities F asymptotically approaches the adi- 
abatic index of a relativistic Fermi gas, T ~ 4/3. These tran- 
sitions occur very near the stellar surface, at r > 0.95r* in our 
evolution coordinates. The dramatic change in F presents a 
kink in the density profile that is difficult to resolve, especially 
with spectral methods. Because the kink occurs over a range 
of radii, it cannot be captured by a subdomain boundary, as 
the well-defined surface of a polytrope can. 

We can often improve the performance of our initial data 
solves by manipulating the EOS minimally. We find we can 



bypass most of the non-smoothness of the stellar profile by 
appending a simple polytrope to the low-density portion of 
the cold EOS. This makes the initial star more amenable to 
spectral representation. Specifically, we use P{p) = Kop^° to 
describe pressures and 



e(p) = 



"° -/»-' + e,Mft 



■r 



to describe energies below pbreak- Pbieak is determined by en- 
forcing continuity in the effective adiabatic index (FCpbreak) = 
Fo), K{) by enforcing continuity in pressure, and eshift by en- 
forcing continuity in energy. We use Fq = 2 with the LS EOS 
at r = 0.01 MeV and Yg in neutrinoless /3-equilibrium. This 
implies pbi-eak = 4.49 X 10'^^gcm"^ ko= 1.40x lO'^g"' cm^s"^, 
and e,Mf. = 9.13 xlO'^ergg"!. 

The effects of this manipulation on the bulk of the star are 
minimal. It is helpful to bear in mind that the central density 
of a 1.4Mo star described by this EOS is ^7 x 10''*gcm"-'; 
we append the polytrope more than an order of magnitude 
below this density. Though the star's radius decreases by sev- 
eral percent, its rest mass changes by one part in 10^. Test 
evolutions of isolated stars show that perturbations from equi- 
librium, due to switching from the modified EOS used in the 
initial data to the true EOS used in the evolution, are smaller 
than the numerical noise at our typical resolutions. 

In practice the manipulated density profile still falls off 
somewhat non-smoothly at the surface, and we can some- 
times gain further accuracy in our initial data by adding an 
additional subdomain, isolating the outer 20% of the star (in 
radius) to a thin spherical shell. 

Note that we limit these two surface-capturing methods to 
the initial data solve, during which the metric and fluid vari- 
ables are represented on spectral subdomains. During the evo- 
lution, the fluid variables are represented on a finite-difference 
grid, where we also employ a high-resolution shock-capturing 
technique that was designed to describe fluid discontinuities. 

Because the quasi-equilibrium approximation ignores the 
infall motion that should be present initially, using this data 
directly results in eccentric orbits, which could affect the disk 
and outflow masses. Therefor e, we perform on e iteration of 
eccentricity reduction (see Pfeiffer et al.||2007 for details) to 
incorporate the initial infall and to reduce the initial eccentric- 
ity to -0.01. 

2.3. Neutrino leakage scheme 

Neutrino leakage schemes attempt to account for the en- 
ergy loss and electron number alteration caused by the emis- 
sion of neutrinos by the nuclear fluid. To accomplish this, 
one first estimates the local effective energy emission rate Q^ 
(ergs"' cm"-'), the effective lepton number emission rate R^ 
(s"' cm"^), and the neutrino radiation pressure P^ as measured 
in a local Lorentz frame comoving with the fluid. In such a 
frame 

d,{pY,) = -R,mu , (2) 

a,r°"=-e, , (3) 

d,Tj,,=-djP, , (4) 

where niu is the atomic mass unit, and T"'^ is the stress tensor. 
Then the energy, momentum, and composition source terms 



follow directly from putting the above in generally covariant 



form. In the comoving Lorentz frame of Equations 2|4 ifi = 1, 



u-i = 0, and gap = r]ap, so a covariant form of the neutrino 
sources is 






-RuU , 



(5) 
(6) 



For the time derivative of P^, we assume simple advection 
u'^V aPv = 0, which is not quite right, since the neutrino field 
changes in ways other than advection, but the effect of the 
diPu source in our simulations (and, indeed, the effect of all 
neutrino pressure terms) is quite small, so a more accurate 
estimate is not necessary. 

For the computation of Q^ and R^, our code is essentially 
an extension of the GRID neutrino leakage code ( |0' Connor] 
|& Ott|[2010| l to three dimensi ons. That code in tu rn c losely 
follows the leakage schemes of |Ruffert et al.| ( |T996| l and |Ross^ 
wog & Liebendorfer (2003]). In optically trim regions, emis- 



sion rates should be given by the local rates of neutrino- 
generating interactions for each neutrino species i/,: 2*'^'' and 
R^l'^''. In optically thick regions, energy and leptons escape via 
diffusion at rates Q**^ and /?f*. The effective rates come from 
interpolating these; e.g., for energy emission. 



free /odiff 



Q.,= 



QTQt] 



diff 



(7) 



Like most other leakage schemes, we include three neutrino 
species, i/^, V^, and i^x, where I^^ includes ji and r neutrinos 
and antineutrinos. Then Q,y = Qu+ Qu +^Qiy and Ri,=Ri, - 

For free emission rates, we include /3-capture processes 
{e~p — > nvg and e^n — > pV,,), electron-positron pair annihi- 
lation, plasmon decay, and nucleon-nucleon Bremsstrahlung 
emission. (For the latter process, we use the rate given in 
jBurrows et al.l]2006] ) In calculating the opacity, we include 
scattering on nucleons and heavy nuclei and absorption on 
nucleons. All particle species, including the neutrinos, are 
assumed to have Fermi-Dirac distribution functions. 

As in GRID, we assume the neutrino chemical potentials 
are 



^,,=Ai:,^(l-e-<^^.>), 



(8) 



where /it^ is the /3-equilibrium value of ^i,. (/i^^ = -p^ = 

fJ-e + f^p- fJ-n, IMa = 0) and (ti,.) is the energy-averaged opti- 
cal depth. 

Given the primitive variables and neutrino chemical poten- 
tials, we can compute t^.{E^.), the optical depth of neutrino 
species j/, at energy E^. along some path parameterized by £, 
with the line integral 



T,,(E,,)= [X,,{E,,)]-'de = El \ 



(9) 



where X^(E^.) = Ei,~^X^~^ is the mean free path, and we have 
taken advantage of the fact that all the absorption and scatter- 
ing cross sections considered have a common neutrino energy 
dependence, a,y. ex El,, which can be factored out. Using the 
optical depth, the diffusion rates Q** and Z?**^ are estimated 
as in Rosswog & Liebendorfer (2003) . 
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Figure 1. Schematic of the optical depth calculation. The fine (dotted) grid 
represents the finite difference grid for the fluid evolution; the coarse (dashed) 
grid represents the opacity grid for the optical depth calculation. The blue 
region in the lower left represents dense fluid in which the neutrino mean free 
path is small. For each zone of the opacity giid, we first integrate five axial 
rays, excluding the —z. axis because of equatorial symmetry; three of these 
are labeled: x™, xP, zP. Then we integrate two additional 'most promising' 
diagonal rays, lying between the minimum optical depth axis and the two 
next-to-minimum axes; one of these is labeled: d'. 

The calculation of t^^ uses interaction cross sections that 
depend on /i^,, and by Equation ]8 these themselves depend 
on (t^,). grid iterates Equation 8] but since optical depth 
calculations are more expensive in 3D, and the interpolation 
formula is itself somewhat arbitrary, we have experimented 

for 



with introducing a second optical depth variable (r, 
use in Equation ]8] We have tried 



I'i I approx 



1) setting 



' Vi j approx 



equal to the energy average of Tj,.(E^.) 



and, as in GRID, iterating to rough convergence 

2) limiting the code to one iteration per opacity calculation, 
starting from previous calculation's results for (Ty,)approx and 
fi,y.. Since the matter distribution changes little between opac- 
ity computations, this method is very similar to the first. 

3) using an analytic fit for neutron stars: 



log 



10\' 1^,7 approx 



= 0.96(logioPcgs-11.7) 



(10) 



4) modifying the above to capture the effect of temperature on 
the average neutrino energy, and hence the cross-section: 



log 



lOV'i'i/ approx 



= 0.96(log,p.,-lL7)(^^) .(11) 



The effect of the assumption for (t^,) approx on the luminos- 
ity is fairly small for all epochs of our simulation (of order 
10% or less), and the runs below use the simplest scheme (3). 
The effect on the lepton number emission rate can be tens of 
percents at times, meaning the choice of neutrino chemical 
potential may have some effect on the composition evolution 
of the disk, an issue we intend to explore more deeply in fu- 
ture work. Numerical experiments do show that the most no- 
table aspect of the disk Yg evolution, the deleptonization at 
late times when the disk cools, is also found with scheme 2 
and does not seem to be an artifact of this chemical potential 
choice. It should be emphasized that (T,y) approx is only used 
for Equation IH] The diffusion rates are determined using t^. 
as computed in Equation ]9] 

The most costly part of the leakage scheme is the estimate 



of the optical depth. To compute this, we first interpolate the 
opacity (factoring out the E^ neutrino energy dependence, as 
discussed above) onto a lower resolution 3D grid. Then we 
compute line integrals of this quantity in 7 directions to the 
computational boundary: one integral forward and backward 
on each coordinate axis, excluding -z because of our grid's 
reflection symmetry, plus two diagonals based on the 'most 
promising' of the coordinate directions (see Fig.fTll. The min- 
imum line integral gives the optical depth. (More precisely, 
the optical depth is E^ times this integral.) 

To test our implementation of the leakage scheme, we have 
constructed spherically symmetric configurations of hot gas 
and compared our results {R,y., Q^., and the optical depth) with 
those of GRID. We set the neutrino chemical potentials in the 
same way for this test to have a clean comparison. We tried 
spheres with a range of central density, temperature, and Y,,, 
so that some were optically thin and others optically thick. As 
expected, the two codes agreed. 

The emission rates also allow us to compute the total neu- 
trino energy and number luminosity radiated to infinity by 
performing proper integrals over Q^ and R,y . For the energy 
luminosity at r — > oo, L^, we multiply the integrand by a red- 
shift factor ^00 : one factor of y^goo for the time dilation and 
one factor for the energy redshift. For the integral of R^, we 
multiply by a redshift factor of y^oo- 

3. SUMMARY OF EVOLUTION: GLOBAL MEASURES AND 
CONVERGENCE 

Our evolution begins with a coordinate separation of 83 km 
between the centers of the two compact objects. The inspi- 
ral takes place over 6 orbits (28 ms) and is followed by tidal 
disruption and merger We continue to evolve the spacetime 
and fluid ^7 ms past merger, by which time the metric is 
mostly stationary, but the accretion disk remains highly dy- 
namic. The interpolation and communication of the fluid vari- 
ables to the spacetime (pseudospectral) grid is the dominant 
computational cost during this phase of the evolution, so con- 
tinuing to evolve the complete system for many tens of mil- 
liseconds would be prohibitive with the current code. Instead, 
we continue the evolution some 40 ms longer in the Cowling 
approximation, that is we evolve the fluid but leave the space- 
time fixed. 

The Cowling approximation ignores several effects. First, 
it neglects changes in the disk's and the black hole's gravi- 
tational pull. The mass accreted onto the black hole is some 
two orders of magnitude smaller than the black hole mass, 
so this effect is probably unimportant. Second, the Cowling 
approximation cannot capture instabilities due to the disk's 
self-gravity. However, this should not qualitatively alter per- 
turbations in the bulk of the disk. The threshold for gravi- 
tational instability can be estimated from Toomre's criterion: 
Qt = ncJiiiGYi) > 1 for stability (Safronov 1960 Toomre 
1964| l, where Cs is the sound speed of the disk, k is the 
epicyclic frequency, S is surface density, and G is the gravi- 
tational constant. For the disk considered here, the minimum 
Toomre parameter there is about 20 except for the inner and 
outer edges. The small amount of matter in the outer regions 
of the grid is not yet in circular orbit equilibrium, so the sta- 
bility condition is inapplicable. At the innermost region of 
the disk, k ^- at the innermost stable circular orbit, but 
here the Cowling approximation does capture the main or- 
bital instabihty. Third, the Cowling approximation discards 



the gravitational waves caused by the disk's nonaxisymmetry. 
For the observed mass quadrupole variations induced in the 
disk, gravitational waves will affect the modes on a timescale 
Emode/Lcw ^ 10^ s (whcrc iimode ^ Mov^((5S/S])^ is a charac- 
teristic energy of the spiral waves, and Lew is the gravitational 
wave luminosity). Thus we may safely ignore this effect. 

To check the accuracy and robustness of our results, we 
evolve the merger at 3 resolutions, which we label 'LI', 
'L2', and 'L3'. The effective number of points on the spec- 
tral/fluid grid are roughly 70^/120^ for LI, 75^140^ for L2, 
and 87''/160^ for L3. Since the spectral evolution uses adap- 
tive meshing, the number of colocation points changes sig- 
nificantly over the course of an evolution, so the above are 
averages. Also, the actual numbers of stored grid fluid points 
are half the above numbers, since each stored point gives the 
fluid variables at a point below and a point above the equator, 
following our assumption of equatorial reflection symmetry. 
After tidal disruption, the coordinate system of our fluid grid 
is driven nonuniformly to higher resolution close to the black 
hole. Our high resolution corresponds to a fluid grid spacing 
of ^160 m during inspiral and about 0.7 km vertical / 2 km 
radial in the inner disk after merger. 

In addition to the convergence tests, we perform one 
merger, at the LI resolution, without neutrino leakage. We 
label this run 'Llno;^'. Comparing Llnoi' to LI lets us iso- 
late the effects of the neutrinos on the fluid. 

In Figures [2] and |3] we show the evolution of several main 
global quantities for all runs. Comparing resolutions allows 
us to discern which evolution features are robust and which 
are numerical artifacts. For example, the neutrino luminosity 
(L^) spikes 2 ms after merger, but this feature appears to con- 
verge away as resolution is increased. Volume renderings of 
the fluid (see an equatorial slice in Fig. |4]i show that at these 
times, some of the nuclear matter is spread in a thin tidal tail 
which is underresolved at low resolutions. The luminosity at 
later times comes from the much better-resolved disk, and we 
see good agreement between the higher resolutions on L^, af- 
ter about 5 ms. Figure |2] reveals that the rest mass of LI at 
late times differs from that of the higher resolutions signif- 
icantly more than the high resolutions differ between them- 
selves. By 10 ms this difference has settled to about 30%. 
The reason for this is that resolution LI generates a hotter, 
fatter, lower-density disk. The puffed-up disk drives more of 
its matter off of our computational domain, an effect that we 
can already see in outflow measurements by 2 ms after merger 
(see Section!?]). Finally, in Fig.[3| we again see that the dif- 
ferent resolutions agree well on the post-merger composition 
and entropy, giving us confidence that these are roughly cor- 
rect, though strict convergence is lacking. Given the strong 
shocks and turbulent-like behavior in the disk, this is unsur- 
prising. However, comparing the leakage runs to Llnoi^, we 
see that the secular neutrino effects are much larger than the 
differences between resolutions. 

4. TIDAL DISRUPTION AND EJECTA 

In a BHNS merger, matter expelled at high velocity may 
ultimately become unbound from the central gravitational po- 
tential. In addition to enriching the inte rstellar medium with 
r-process elements (iL attimer & Schram m|1974 Freiburghaus 
|etaL|[T999i fAmould et al.||2007| [Korobkin et al.||2012| i, this 
nuclear matter may emit an electromagnetic signal. The ra- 
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Figure 2. Global convergence measures, showing the three resolutions with leakage, (merge is the time at which 50% of the rest mass has depleted. We use 
different scalings for the time-axis before and after 8 m s in order to highlight the initial disruption and accretion. Top panel: L- norm of global generalized 
harmonic constraint violations (seelLindblom et al. 12006). We stop tracking this measure at ~7 ms when we freeze the spacetime metric and begin the Cowling 
evolution. Middle panel: baryonic rest mass remaining on the computational grid (in units of Mq). The initial decrease beginning at —1 ms is driven by accretion 
onto the black hole. The second decrease beginning near 3 ms is due to matter falling off the outer bondary of the computational domain. The disk mass quoted 
throughout (Mq Ri 0.3Mq) is the gravitationally bound mass outside the BH at 5 ms in the L3 evolution, i.e. the mass on the grid minus the instantaneous estimate 
of the remaining unbound mass discussed in SectionWl Bottom panel: total neutrino luminosity (in units of 10^' erg s"' ), as measured at r — > oo. 
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Figure 3. Secular chemical and thermal evolution of the disk for the three resolutions with leakage, and also hlnou. Top panel: mass-weighted average of the 
electron fraction. In the Llnoi/ evolution, Yg simply advects with the flow. The early decreases at ms and 3 ms are due to the major mass loss episodes: initial 
accretion, and loss of the tidal tail from the outer boundary. Bottom panel: mass-weighted average of the entropy, in units of kg baryon"' . 

amplify and accelerate electrons and positrons, thus emitting 
synchrotron radiation. This radio emission continues for a 
deceleration timescale (months to years) dependent upon the 
density of the surrounding medium (|Metzger & Berger|2012{ 
INakar & Piran|201 l[|Piran"et al.|2013| l. 

These transient electromagnetic signals are of sufficient in- 
terest to warrant a thorough examination of any ejecta in our 
simulation. There are two dominant conveyers of ejecta in 
BHNS mergers: tidal tails and accretion disk winds. In this 
simulation, most of the unbound mass is produced by the 
tidal tail. A smaller amount is ejected in a plume during 
the early merger phase, as the infalling gas stream collides 



dioactive nuclei formed in the neutron-rich fluid quickly fis- 
sion, emitting high-energy beta- and gamma-radiation. This 
heats the ejecta, which becomes optically thin on an expan- 
sion timescale (days to weeks), and radiates thermalized pho- 
tons, producin g an isotropic kilonova (or ' macronova', see Li 
& Paczynski|[T99 8 ; Roberts et al. 201 \\ [Metzger & Berger 
2U12[ | Rosswogl[2012nPiran et al.||20l3t (Kasen et al ||2013| l 
at optical, or perhaps infrared, wavelengths, depend ing on 
the highly-u ncertain opacity properties of the material ( [Kasen 
|etal.|2013| l. 

Additionally, the high-velocity ejecta form a blast wave in 
the circumbinary medium. At the shock front, magnetic fields 
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and shocks with itself. Also a small amount of outflow from 
the late disk is observed, but at low densities where numeri- 
cal errors can introduce spurious acceleration and heating of 
the fluid. Note that this simulation ignores some effects that 
could drive strong winds (e.g. neutrino heating) and does not 
evolve long enough to see some others (e.g. large-scale He re- 
combination). Below we review previous results, describe the 
formation of the tidal tail, describe our method of analyzing 
the ejecta, and finally report measures of its mass (Mgj), ki- 
netic energy (fiq), distribution of velocities, and composition 

((i;)ej). 

There are a number of recent s tudies examining ejecta fro m 
binary neutron star mergers (e.g. Rosswog 2012; Hotokezaka 



et al.| 201 3i and high-eccen t ricity mergers (e.g. Ste phens, 
et al.poill [East et al.||2"0T2i [East & Pretor ius 2012). But 



here we survey results from low-eccentricity BHNS stud- 
ies to emphasize t he infl uence of black hole spin. |Lat-| 
[timer & Schramrn] ( |1974] ) made semianalytic estimates in a 
Schwarzschild spacetime. They found Mq ^ to O.14M0 for 
stars that disrupt close to the blac k hole; no mass is ejected 
from stars that disrupt far away. |Rosswog ( 2012[ l simulated 
two cases of ^ ^ 4 and ^ ~' 7 in a Newtonian potential. He 
found Mej < 0.05 M0. Relativistic simulations have charac- 
teristically yielded more conservative ejecta estimates. No- 
tably, high mass-ratio, compact neutro n star, low black hole 
spi n systems do n ot even disrupt. (See ,Pannarale et al.|2011 
and Foucart 2012 f or phenomenological m odels covering this 
parameter space.) Kyutoku et al. ( 2011| l simulated a suite 
of tens of mergers with mass ratios of q ~ 2 to 5 and BH 
spins of a* = -0.5 to +0.75. They noted the possibility of 
Mej ~ 0.01 M0 for systems with prograde BH spin and stiff 
EOSs. Recently, Foucart et al. (2013 ) simulated several merg- 
ers of q = l, with large BH spins of \a* \ = 0.9 varying in incli- 
nation. In cases of aligned spin and orbital angular momen- 
tum, they found ejected masses of Mej ~ 0.09Mq. Finally, in 
a study with nearly -extremal BH spin, [Lovelace et aLl ( |2013| l 
found Mej could be as high as 0.3 Mq. It appears, then, that 
for BHNS systems, the region of interest for significant ejecta 
may be the parameter region of high spin. 

In the present simulation, we see a large tidal tail form just 
before merger. About two orbits before merger (t = -3 ms) 
the coordinate separation of the two centers of mass has de- 
cayed to 40 km and the neutron star has become extremely 
distorted. After another revolution (t = -1.5 ms) the separa- 
tion has decayed to 25 km, the star overflows its Roche lobe, 
and it begins to accrete onto the black hole. A tidal tail has 
already formed from the trailing edge of the star, extending 
outward and lagging an orbit behind the core. At f = 0, the 
core falls into the black hole. Over the next 3 ms the tidal tail 
sweeps out and away from the black hole (see Fig. |4]i. We 
follow the tail by periodically resizing our computational do- 
main outward to a radius of ^400 km, where we allow the 
tail to fall off of the grid. From its formation until its exit, the 
tidal tail is evolved on the computational domain for ^5 ms. 

For a stationary spacetime, m, (the projection of the 4- 
velocity along the timelike Killing vector field) is a constant 
of motion for geodesies. Assuming it is also asymptotically 
flat, u, gives the asymptotic Lorentz factor (7 = -u,) of the 
geodesic as it escapes to infinity. To the approximation that 
the spacetime is settled and that pressure is dynamically neg- 
ligible, u, will be a constant along fluid parcels, and matter 
with M, < -1 may be flagged as unbound. Neither approxi- 



mation is strictly true, but both become better satisfied as the 
outflow expands outward away from the dynamical region and 
decompresses to low pressures. Thus, u, of fluid elements in 
the tail, and the total amount of unbound material in the tail, 
should become constant as the tail expands. If this settling 
happens before the tail leaves the computational grid, mean- 
ingful statements about the amount of outflow and its asymp- 
totic velocity distribution can be made. We integrate the mass 
satisfying the unboundedness condition u, < -I over times 
from f ^ ms, just as the tidal tail is forming, to f ^ 7 ms, 
after the tail has entirely fallen off the computational domain. 
We mask out regions within 5OM0 (73 km) of the BH in the 
evolution coordinates. At this distance the Newtonian gravi- 
tational potential is Msn/r ;< 0.15. The integral of total mass 
ejected is robust against 30% variations in the mask radius. 

Measures of Mej and E^j from the three resolutions are 
shown in Fig. [5] where we have calculated the total integral 
over unbound matter on the grid. This raw measure increases 
early as a growing amount of material is flung out from the 
disk and exits the masked region. In addition to the tail, vol- 
ume renderings of the unbound matter reveal the existence of 
a plume of outgoing matter ejected above the equatorial plane, 
of lower density than the tidal tail visible in Fig. [4] This sec- 
ond outflow is produced during the merger, pemaps by the 
shock waves in the infalling matter. The most energetic mat- 
ter in this plume begins to fall off the computational domain 
after 2.5 ms, whereas the dense tidal tail begins to fall off after 
3 ms. Therefore, we use the peak values of mass and kinetic 
energy in Fig.|5]as conservative estimates. We use Richardson 
extrapolation to derive uncertainties from the highest two res- 
olutions, and we assume 2"'' order convergence, though Fig.pl 
indicates approximately 6* order covergence. This method 
gives Mej = 0.08 ± 0.07 Mq, and ^ej = 10± 17 x 10^' erg. 
These errors are upper bounds which are probably large over- 
estimates. Nonetheless, even at our highest resolution, the tail 
is poorly resolved: at ~'70 km wide, it is covered by ~'15 
grid cells. This is a common limit of grid-based methods; 
smoothed-particle hydrodyanmics methods neatly overcome 
this limit, but find their handicap in evolving accurate initial 
conditions for the tail. 

The velocity distribution peaks across all resolutions at 
0.2c, but the dispersion in velocities is large, with some mat- 
ter escaping at v > 0.5c (see Fig. |6|. Because average tem- 
peratures in the tidal tail are ^1 MeV as it leaves the com- 
putational domain, the fluid is still in nuclear statistical equi- 
librium, so our equation of state is still valid. The material 
ejected from the system is neutron-rich, with average electron 
fraction {Ye)ej «0.1. 

5. RELEVANT TIMESCALES FOR THE DISK 

The density in the accretion disk peaks at a circumferential 
radius r of around 50 km from the black hole, and most of the 
disk mass is within 100 km of the black hole. At the density 
peak, the density is p « lO'^g cm""*, the angular frequency ft 
is around 3 kHz, and the sound speed c, is around 0.1c. The 
disk is hot and thick: H/r ^ cJiVlr) ^ i. The dynamical 
timescale of the disk is, to order of magnitude. 
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Angular frequency decreases with distance from the black 
hole, so the orbital period at r =120 km (roughly speaking. 
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Figure 4. Equatorial slices of rest density from L3 during disruption and tidal tail formation at 1, 2, and 3 ms after merger (from left to right). The unbound 
material is outlined in black. The field of view is fixed in all panels, so that each has a scale of approximately 680 X 680 km^. The computational grid expands 
with the tail until it covers a circle of radius ~400 km (in our asymptotically flat evolution coordinates) at which we remove matter from the grid. Note that we 
distort the coordinate ma p between the fu ndamental evolution frame and the frame in which the finite difference grid is uniform in order to concentrate resolution 
near the black hole (see |Duez et al.|2010) . 

' ' ' ' ' ' ' the disk's outer edge) is 10 ms. Thus, even this outer region is 

evolved for several dynamical times. The perturbations in the 
inner disk rotate at frequencies similar to the local flow and 
thus orbit the black hole on millisecond periods. 

The disk has strong entropy gradients in the inner and outer 
regions, so disk gas can experience buoyancy forces there on 
a timescale given by the Brunt- Vaisala frequency A^bv, 
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Figure 5. Integrals of the unbound fluid on the grid. Top panel: total asymp- 
totic kinetic energy. Bottom panel: total rest mass. For both integrals a sphere 
of 73 km near the black hole was masked out. The drop-off after 3 ms is due 
to the loss of matter through the outer boundary. 
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where g = p ^ VP and S is the specific entropy. 

In addition to these dynamical processes, there are secular 
processes that alter the disk structure on longer timescales. 
One such effect is the depletion of baryonic mass in the disk. 
The disk's mass is Mq ~ O.3M0 after merger (see Fig.l2]i. For 
the first 30 ms after merger, accretion onto the black hole pro- 
ceeds at a rapid rate of Mq w 2M0 s~'. During this time, the 
gas in the outer disk settles into circular orbit, which requires a 
transfer of angular momentum from the inner disk. (See Sec- 
tionl6]for a fuller discussion of the disk's angular momentum 
evolution.) The transport of angular momentum away from 
the inner disk by spiral density waves drives accretion onto 
the black hole. Eventually, the high-density middle region 
of the disk settles to equilibrium, but dynamical, nonaxisym- 
metric distortions persist near the disk's inner and outer edges 
(also discussed in Sectionl6]l. Spiral waves thus travel outward 
and drive a reduced rate of accretion throughout the simula- 
tion. After 30 ms, the accretion rate has dropped by nearly an 
order of magnitude. At these late times, disk mass is also lost 
to a weak outflow from the inner disk at a comparable rate to 
the accretion into the black hole. Combining these effects, we 
can define a mass depletion timescale 
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Figure 6. Unbound rest mass binned by magnitude of the asymptotic fluid 
velocity. We show LI, L2, and L3 at ~2.7 ms, just before the peak in Fig. [5] 



The disk's thermal evolution is driven by shock heating, 
compression, advection of heat into the black hole, and ra- 
diative cooling. The disk's thermal energy is defined as 



P*[^-^co\i{P,ye)]d X, 



(15) 



where e and ecoid(p,J'(>) are the actual specific internal en- 
ergy and the specific internal energy at the lowest temper- 
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ature in the EOS table, for which the gas is degenerate. 
10^^ erg. The total neutrino luminosity is 
ergs"'. Therefore, the cooling timescale due to 



We find ^thermal 

u ^ io"-54 

neutrinos is 



logioQv (erg Ion's"') 
-47 



^cool ' 



^thermal 



lOms-lOOms. 



(16) 



It is important to remember that the disk is partly pressure- 
supported. Thus radiative energy loss may come from the po- 
tential energy reservoir via disk contraction, in addition to the 
more obvious thermal energy reservoir via disk cooling. In- 
deed we find that evolving without neutrino cooling leads to a 
disk that is not only hotter but also much more extended and 
less dense. The timescale for composition change is Ne/Rv, 
where Ne is the number of electrons in the disk and R^ is the 
total net lepton number change rate due to neutrinos. This 
is initially also about 10 ms. Then, 20 ms after merger, bal- 
ance between v,, and Ve emission is roughly achieved, and 
the composition subsequently changes more slowly. Thus, 
the neutrino emission significantly influences the energy and 
composition of the disk over its lifetime. 

6. DISK FORMATION AND FALLBACK 

As can be seen from Fig. [3] as the accretion stream collides 
with itself, shocks heat the gas for roughly one millisecond, 
until the density-averaged entropy settles at ~8A;Bbaryon"'. 
A hot accretion disk forms in the vicinity of the black hole. In 
Fig. 17] we show density snapshots of the disk at a represen- 
tative time ^30 ms later, and in Fig. [8] we show azimuthally 
averaged equatorial density as a function of circumferential 
radius. The maximum density remains at a fairly steady level 
of ^lO'^gcm""*. As shown in the bottom panel of Fig. ^ the 
densities and temperatures are sufficient to render the interior 
of the torus optically thick to all species of neutrinos, with 
optical depths (averaged over neutrino energy) of order 10. 

The disk is initially extremely distorted and nonaxisymmet- 
ric. For a completely stable disk, one would expect the inner 
regions, where the dynamical timescale is shortest, to settle to 
a stationary axisymmetric state before the outer regions — as 
was seen, for example, in a recen t BHNS F-law EOS m erger 
carried out with the same code ( [Lovelace et al.| |2013). We 
see instead that it is only the middle disk that settles to an ap- 
proximately axisymmetric state after about 30 ms. In the inner 
disk, very close to the black hole, clear, order unity deviations 
from axisymmetry in the form of trailing one-armed spirals 
persist throughout the evolution (see the bottom panel in Fig. 
17]). These perturbations appear at the disk's inner edge, rotate 
at roughly the rate of the local fluid (^ms periods), expand 
outward and dissipate, and then reform many times during the 
disk evolution. Such behavior suggests that a fluid instability 
may be preventing the disk from promptly settling. 

Further insight into the structure of the disk comes from the 
profiles of the specific orbital energy {E = -u, - 1) and angu- 
lar momentum (L = -u^/u,) shown in Fig. ]9] Each panel in- 
cludes three curves. One is the actual E or L profile, measured 
on the equator and averaged over azimuthal angle. Another is 
the 'Keplerian' E and L, the values that would be found for 
geodesic equatorial circular orbit given the evolved spacetime 
metric of the system. The geodesic circular orbits are those 
for which WffU = 0. These curves thus include the effects of 
the disk's self-gravity (at least, as it was at the beginning of 
the Cowling evolution) but disregard pressure forces. We see 




Figure 7. Three-dimensional distribution of neutrino energy loss and fluid 
density from L3, 30 ms after merger, depicted in evolution coordinates at a 
slanted view. Top panel: effective local neutrino power, Qi, (no redshift ap- 
plied), summed over all species (in units of ergkm"^ s"' ). Middle and bottom 
panels: meridional and equatorial slices of density in the fluid rest frame. 
The equatorial slice reveals a spiral mode. In all three panels, densities below 
lO'^gcm"^ have been masked out to show the structure of the disk. The disk 
radius and half-thickness are 1 10 km and 45 km, respectively. 
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Figure 8. Top panel: azimuthally-averaged profiles of the density p togetiier 
witli its ]?MS deviation from axisymmetry, 5p,n,s = ((P~(p))~)''^> plotted as 
functions of tiie circumferential radius at the equator r. Values for 17 ms and 
31ms after merger are shown. Bottom panel: density- weighted azimuthal 
and vertical average of the energy-averaged neutrino optical depth for each 
species of neutrino, computed 3 1 ms after merger 

that the angular momentum profile, for instance, is much shal- 
lower than would be found for an equilibrium thin (geodesic) 
disk. Finally, we include 'equilibrium' curves, the E and L for 
equilibrium circular orbit given the existing metric, density, 
and pressure. These are orbits for which WffU = -WP/{ph). In 
regions where deviations from equilibrium are nonlinear, one 
must take care in identifying these curves with the true equi- 
librium about which the disk is perturbed, since m = per- 
turbations in the pressure will feed back into the equilibrium 
condition. Given this qualification, we see that the highest- 
density regions do appear to be in equilibrium in an angle- 
averaged sense. The outer disk has sub-equilibrium angular 
momenta, so the gas remains in eccentric orbit. Interestingly, 
the energy and angular momentum increase dramatically at 
the inner edge of the disk. This feature is not present in the 
geodesic curves, but it is in the equilibrium curves, indicating 
that it is a consequence of a sharp pressure gradient near the 
inner edge of the disk. Below, we will consider the effect of 
these gradients on the expected stability of the disk. The high 
energies in the inner disk would make it easy to generate an 
outflow there. We do indeed observe some mass ejection from 
the inner disk, but at densities too low to be reliably modeled 
by our code. 

The difference between the actual and equilibrium energy 
curves provides a measure of the mode energy, again, to the 
extent that equilibrium curves can reliably be computed from 
an incompletely settled pressure profile: Smode = J P*{E - 
£^equiiibriuni)<^^^-'f- The modc energy is positive but decreases 
during the early, rapid accretion, phase. Then it settles and 
even shows episodes of growth, which may indicate that in- 
stabilities are stimulating these modes. 



In Fig. 10 we plot the density-averaged pressure P, entropy 
S, temperature T, and electron fraction F^ as functions of cir- 
cumferential radius — reducing from three dimensions by av- 
eraging these functions vertically and azimuthally. The en- 
tropy profile is fairly flat in the bulk of the disk, with the ex- 
ception of a steep negative gradient (dS/dr < 0) in the hot 
inner region. Since in this region pressure increases with ra- 
dius, this entropy gradient is actually a stabilizing force. The 
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Figure 9. Energy and angular momentum of the disk 3 1 ms after merger. Top 
panel: azimuthally-averaged equatorial orbital energy E = -u, — 1 . Bottom 
panel: azimuthally-averaged angular momentum L = —u^/ut. In addition to 
E and L for the actual flow, we include values for circular orbit geodesic 
motion and circular orbit equilibrium (i.e. including pressure) motion. 




r (km) 

Figure 10. Profiles of the azimuthally- and vertically-averaged pressure, 
entropy, temperature, and electron fraction 3 1 ms after merger. For the pres- 
sure, we plot separately the contributions from the neutrinos and everything 
else (nucleons, electrons, photons). For {Y,,), we show both the actual value 
and the value for neutrinoless /3-equilibrium for the given density and tem- 
perature profiles. 
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Ye gradient also affects the buoyancy, but for these profiles, its 
effect is much smaller than that of the entropy. The relativistic 
Solberg-H0iland critera for convective stability of an axisym- 
metric equilibrium fluid are given, f or a coordinate basis in 
which g„. = g,.^ = 0, by ( |Seguin|1975[ ) 



10 



7 • VL+(phy\dU/dS)pVP • V5 > , 
(dU/dS)p(^ X VP) • (V5 X VL) > , 
where U = p+peis the total energy density and 



7 = (m''mo)2[(1- 



^UQ^WL-vn] 



(17) 
(18) 

(19) 



Equations 17 and fTSl correspond to the familiar Newtonian 
radial and vertical stability conditions, respectively. The first 
term in Eq. 17 dominates everywhere inside of r w 120 km, so 
dL/dr determines the stability to radial convection through 
most of the disk. The condition for stability is met except 
in the inner edge (r < 27 km) of the disk. This inner region 
should be unstable. The unstable gradient persists, and it is 
presumably connected to the inability of the inner disk to set- 
tle to stationary axisymmetry. The outer disk also has an un- 
stable angular momentum gradient (dL/dr < 0), but it is stabi- 
lized by a strong positive entropy gradient. The entropy does 
not change much with height, but we do identify a negative 
gradient in a region near the equator, leading to an instability 
according to Eq. 18 there. Buoyant forces would be expected 
to correct this gradient on timescales (-N^y z)"'''^ ~ 10ms, but 
apparently other processes (see SectionlTlon processes that af- 
fect entropy evolution) are sufficient to keep the entropy pro- 
file roughly fixed. 

The temperature is high (^-^8 MeV) in the inner disk, and 
in these regions, the contribution of trapped neutrinos to the 
total pressure is as high as 12%. In the bulk of the disk, the 
neutrino pressure is negligible. 

7. NEUTRINO-DRIVEN DISK EVOLUTION 

From Fig. |2] we see that the late-time neutrino luminosity 
is around 2 x 10^^ ergs"' and continuing to drop. During the 
early disk phase, M « IM0S"', so the accretion efficiency 
is ?7 = L^/Mc^ « 0.1. The disk remains very nondegenerate 
throughout the simulation, with the thermal component of the 
internal energy remaining nearly a constant 85% of the total 
internal energy throughout. The actual values of the thermal 
and internal energies decrease at a rate only about half L,y (see 
Fig. \TT\ . The rate of loss of thermal energy due to flow out 
of the grid, comprised at early times primarily of accretion 
into the black hole and at late times primarily of disk outflow, 
is also comparable to L^. Thus, these energy sinks are be- 
ing countered by heat sources sufficient to slow the cooling 
rate by about a factor of 3. Two physical sources of heating 
in parts of the disk are adiabatic compression (primarily the 
observed flattening of the disk toward the equator as it cools) 
and shock heating, the latter presumably occuring in the non- 
linearly perturbed regions of the disk. Adiabatic flattening 
alone is insufficient, since we observe that the rate of total en- 
tropy decrease is also much lower than the radiative entropy 
loss rate. Shock dissipation is the more plausible heat source, 
since losses in spiral mode energy are of the needed magni- 
tude. A third physical energy source, recombination of nu- 
clei, does not occur during the phase of evolution simulated; 
the composition of the disk remains > 99% free nucleons. 
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Figure 11. The thermal evolution of the disk during its cooling phase. L^ is 
the neutrino luminosity, ((^fithcrmal/^Ooutflow the heat loss due to flows out of 
the inner and outer boundaries, and (rf£(hcrmal/^f)disk the numerical derivative 
of the total thermal energy in the disk. (The total internal energy is almost ex- 
actly I.15£theimal thi'oughout the Cowling evolution.) Numerical derivatives 
and outflow measures are noisy, so the curves have been smoothed in time by 
convolving with a Gaussian of width 1 ms. 

A nonphysical source of heating would be numerical viscos- 
ity. This can only be distinguished from true shock heating 
by comparing heating rates at different resolutions. We find 
that L2 and L3 have very similar thermal histories, while LI 
cools 30% more slowly, largely due to a difference in the out- 
flows, although effects of stronger numerical energy dissipa- 
tion might also contribute. The average specific entropy of the 
disk provides a clearer sense of whether the remaining matter 
is heating or cooling. Its evolution is included in Fig. |3] We 
see that, after the initial strong shock heating at disk forma- 
tion, the entropy begins a slow decrease. The fact that differ- 
ent resolutions show good agreement on the entropy evolution 
reassures us that numerical dissipation is not dominating the 
thermal evolution of the disk. 

The disk evolved without neutrino cooling develops a sig- 
nificantly higher, and continually growing, entropy. This 
leads to a larger, more diffuse accretion disk; the maximum 
density is a factor of ten lower than for the neutrino-cooled 
disk, and the density profile lacks the sharp peak in the in- 
ner region seen in Fig. |8] The larger extent of the non-cooled 
disk causes more matter to escape through the outer bound- 
ary. During the first 10 ms after merger, the accretion rate 
into the black hole for cooled and non-cooled disks is compa- 
rable, and so the Lino;/ disk comes to have a baryonic mass 
lower by about 40% than the LI disk. Again, this is a con- 
sequence rather than the cause of the lower density. At later 
times, the accretion rate of the Llnoi/ disk drops to very low 
values, but outflow continues. 

Because of its lower density, the non-cooled disk actually 
has a lower average temperature than the cooled disk (4 MeV 
vs. 5 MeV), despite its higher entropy, during the first 45 ms 
after merger. After this, the temperature of the LI disk drops 
below that of the Llnoi^ disk because of neutrino cooling. As 
expected, neutrino emission strongly affects the evolution of 
Ye. For the Llnoi^ disk, the electron fraction can only change 
because of the accreted or outflowing matter having Y,, that 
differs from the average, and this turns out to cause only small 
changes to the average (see Fig. [3]). 

The Llnoi^ disk shows the same persistent perturbations 
from axisymmetry that are seen in the neutrino-cooled disks. 
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Figure 12. Neutrino emission cliaracteristics by species. Top panel: lep- 
ton number luminosity, R^. (in units of 10^' s"'). Middle panel: energy 
luminosity, L^. (in units of lO'^ergs"'). The total luminosity is given by 
Li/i, +Li7^ +ALu^. Bottom panel: average neutrino energy (£1/,) (in units of 
MeV). All plots are taken from the L3 simulation, and are calculated for an 
observer at r — ^ 00. 

High-energy, non-equilibrium flow in the inner disk is still 
seen. This confirms that cooling is not the main driver of per- 
turbations in the cooled disks. We find that the geodesic and 
equilibrium angular momentum curves show greater devia- 
tion for the Llnoi/ disk — the geodesic curves are essentially 
the same, but dL/dr is about 25% shallower for the equilib- 
rium curve in the bulk of the non-cooled disk — indicating that 
pressure support makes a stronger contribution to the equilib- 
rium of this disk. The entropy profiles have the same general 
shape, and the expected stability properties are similar. 

Fig. [12] breaks down the energy and number emission by 
neutrino species. Initially, Vf. dominates over v^ emission, 
as would be expected for a neutron-rich, proton-poor gas. 
As emission continues, protons become more numerous, and 
the neutrinosphere cools, so that positrons become less com- 
mon than electrons. (The electrons are mildly degenerate, 
l^e/ksT « 1, SO the relative ratio of electron to positron den- 
sity is quite sensitive to temperature.) Thus, R^^ and /?p, be- 
come closer. At 20 ms after merger, these emission rates are 
sufficiently balanced that Y^ thereafter evolves on the slightly 
longer timescale on which the disk itself is changing. In 
Fig.[TOJ we compare the actual F^. to the value for neutrinoless 
/3-equilibrium (Je = Yg^ such that //^^ = = He + fip- iJ„). At 
early times, Yg < Yg^ and Yg is growing. At late times (as in the 
figure), Yg > Yg^ and Yg is decreasing. The outer disk radiates, 
and thus responds to emission rate imbalances, much more 
slowly. The average energy per neutrino, given by L^./R^^, 
averaged in time over the disk evolution, is about 12, 15, and 
19 MeV for Vg, Vg, and i^x neutrinos respectively; the average 
neutrino energies are not constant, but decrease at a rate of 
1 MeV per 10 ms. When one adds together all four species 
of i/x neutrinos, there are still fewer of them emitted than i/g 



or Vg neutrinos, but their average energy is sufficiently higher 
(i/f are emitted from the hotter regions in the disk interior) that 
their combined luminosity is slightly larger than L^^ and L^,. 

8. CONCLUSIONS AND FUTURE WORK 

In this paper, we present simulations of a black hole- 
neutron star binary merger event using a hot nuclear theory- 
based equation of state with neutrino cooling. This allows us 
to address features of the merger — composition changes, neu- 
trino energy extraction — that are inaccessible to models with 
any less degree of realism. While important simulations with 
hot EOS and neutrino leakage have been done in the past, we 
perform them for the first time in full general relativity. We 
are thus able to include distinctly relativistic features, such as 
a high black hole spin. 

For this high-spin, moderate-mass system, we find that a 
significant amount of matter evades prompt capture by the 
black hole. Roughly 0.08 Mq of low-F^ tidal tail material ap- 
pears to be dynamically ejected from the system at mildly rel- 
ativistic speeds. This material will contribute to r-process en- 
richment of the interstellar medium, and may produce quasi- 
isotropic electromagnetic signals in the form of an optical or 
infrared kilonova and a radio afterglow, although our simu- 
lations do not model such processes. Because the opacity 
properties of the exotic nuclei formed in the tail are so poorly 
constrained (although see Kas en et al.| |2013 for progress on 
this problem), quantitative measures"of kilonova detectabihty 



are subject to large uncertainties. Metzger & Berger 
do, however, define a quantitative figure of merit for tl 
tectability of the radio transient. 



FOM, 



/ ^e' 



rad ■ 



\ 10^" erg/ V Icm 



7/8 



11/4 



(20) 



where n is the number density of the circumbinary medium. 
For the ejecta measured in this simulation, and assuming an 
optimistic circumburst density of 0. 1 cm"^ (Berger et al. 2005] 
Soderberg et al.|2006 !), we roughly estimate FOMrad « 0.15, 
which approaches the detection threshold of 0.2 given for the 
Expanded Very Large Array (EVLAyl at 1 GHz, for events 
within a range of 200 Mpc. 

Of the matter that remains bound, 0.3 Mq forms a hot, 
neutrino-optically-thick accretion disk. We find that neutrino 
cooling has a significant effect on the disk structure; uncooled 
disks are significantly more extended. Nevertheless, the disk 
remains hot and thick, with markedly non-Keplerian orbital 
profiles. The innermost disk configuration appears to be un- 
stable, and order unity deviations from equilibrium and ax- 
isymmetry, in the form of millisecond-period spiral features, 
persist for tens of milliseconds after merger. The neutrino lu- 
minosity after merger is initially very high (10^"* ergs"'), but 
disk depletion due to accretion and outflow, and energy de- 
pletion due to radiation cause the disk to quickly dim. Its 
luminosity is able to remain above 2 x lO^-'ergs"' for about 
50 ms. 

Assuming energy deposition efficiencies o f 0.2%-0.5% 
suggested b y prior studies of vV annihilation (|Ruffert et al. 
1997; Birkl et al.||2007l [Dessart et al1[2009l [Harikae et aL 
2010), this simulation would suggest difficulties in powering 
a GRB for timescales much longer than 50 ms, given similar 
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initial conditions. Even so, the GRB signal itself might be 
longer lived, since velocity differences b etween the front and 
tail ends can radially stretch the outflow ( Aloy et al. 2004[ ). In 
fact, worries about the signal lifetime are not yet warranted, 
since the drop in luminosity is most likely an artifact of the 
absence in this simulation of physical processes expected to 
dominate the long-term evolution of the disk. Most obviously, 
accretion disks like thi s one will be subject to the magne- 
torotational instability ( [Balbus & Hawley|[T99T] l. The MRI 
can drive vigorous turbulence which, through viscosity or re- 
connection at small scales, can reheat the disk. Independent 
of its possible effects on the thermal power available, mag- 
netic forces could make possible Poynting flux-dominated 
outflows, tapping either the rotational energy of the disk or 
the spin energy of the black hole. 

Although it is undoubtedly interesting, the system studied 
here is only one case of BHNS merger, and perhaps not a typ- 
ical one. It will be important to carry out simulations of more 
systems, covering a range of neutron star masses, black hole 
masses, and black hole spins. We believe our code captures 
enough of the relevant physics to be suitable for such explo- 
rations; it should certainly improve upon the poly tropic neu- 
tron star models we and others have previously used for such 
surveys. Simulations with this level of realism could also be 
used to study the effects of the different choices of hot EOS 
used to simulate a merger Several hot nuclear EOSs have 
been made publicly available and are suitable for codes like 
oui-s ^ Lattimer & Swesty'1991','Shen et al.'l QQSllShen et al.| 
[20TTp empel et al. 2012; Steiner et al. 2012i 

We hope to make several improvements to future simula- 
tions. Much could be gained for nucleosynthesis estimates 
by tracking the tidal tail out to farther distances and with 
higher accuracy. To simultaneously resolve outgoing ejecta 
and the incipient accretion disk is a challenging task for any 
grid-based code. The solution will likely involve some com- 
bination of increased computational resources and improved 
modeling techniques, such as adaptive meshing or Lagrangian 
(SPH) evolutions of the outflow. We are also working to in- 
corporate missing crucial pieces of physics. Magnetic effects 
require magnetohydrodynamic simulations, while the effects 
of neutrino heating will only be accessible to our models when 
we replace the leakage scheme with some form of genuine ra- 
diation transport. 
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